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We present long-term (~ lO'^M) axisymmetric simulations of differentially rotating, magnetized 
neutron stars in the slow-rotation, weak magnetic field limit using a perturbative metric evolution 
technique. Although this approach yields results comparable to those obtained via nonperturbative 
(BSSN) evolution techniques, simulations performed with the perturbative metric solver require 
about f/4 the computational resources at a given resolution. This computational efficiency enables 
... us to observe and analyze the effects of magnetic braking and the magnetorotational instability 

' (MRl) at very high resolution. Our simulations demonstrate that (f) MRI is not observed unless 

, the fastest-growing mode wavelength is resolved by > 10 gridpoints; (2) as resolution is improved, 

. the MRl growth rate converges, but due to the small-scale turbulent nature of MRl, the maximum 

growth amplitude increases, but does not exhibit convergence, even at the highest resolution; and 
(3) independent of resolution, magnetic braking drives the star toward uniform rotation as energy 
^ ■ is sapped from differential rotation by winding magnetic fields. 
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■ In differentially rotating neutron stars, an initially weak magnetic field will be amplified by processes such as 
0^ magnetic braking and the magnetorotational instability (MRI) [l|, |2| , causing a redistribution of angular momentum. 

Such differentially rotating stars may arise from the merger of binary neutron stars |3j Uj or from collapse of 
massive stellar cores, even if the cores spin uniformly at the outset ^] (see also 0). 
, To better understand how magnetic braking affects differentially rotating configurations, Shapiro performed a 
• purely Newtonian, magnetohydrodynamic (MHD) calculation in which the star is idealized as a differentially 
rotating, infinite cylinder of homogeneous, incompressible, perfectly conducting gas (see also 0). The magnetic field 
Q is taken to be radial initially and is allowed to evolve according to the ideal MHD (flux- freezing) equations. This 
calculation demonstrates that differential rotation generates a toroidal magnetic field, which reacts back on the fluid 
c/3 flow. Without viscous dissipation, the toroidal fleld energy and rotational kinetic energy in differential motion undergo 
periodic exchange and oscillations on the Alfven timescale. The magnitude of these oscillations, and the maximum 
^ , field strength, are independent of the initial magnetic field strength; only the growth and oscillation timescale depend 
on the magnitude of the seed field. If viscosity is present, or if some of the Alfven waves are allowed to propagate out 
of the star and into an ambient plasma atmosphere, the oscillations are damped, driving the star to uniform rotation. 

Cook, Shapiro, and Stephens (IQ] later generalized Shapiro's calculations for compressible stars. In their model, 
the star is idealized as a differentially rotating, infinite cylinder supported by a polytropic equation of state. They 
performed Newtonian MHD simulations for differentially rotating stars with various polytropic indices and different 
initial values of r/|VK|, where T is the rotational kinetic energy and W is the gravitational potential energy. They 
found that when r/|VF| is below the upper (mass-shedding) hmit for uniform rotation, /3max, magnetic braking results 
in oscillations of the induced toroidal fields and angular velocities, and the star pulsates stably. However, when T/|T/F| 
exceeds /3max, their calculations suggest that the core contracts significantly and shock waves are generated, while the 
outer layers are ejected to large radii to form a wind or an ambient disk. 

Liu and Shapiro [ij carried out both Newtonian and general relativistic MHD (GRMHD) simulations on slowly 
and differentially rotating incompressible stars. They considered the situation in which T <^ A4 ^ \ W\, where Ai is 
the magnetic energy. Due to the assumptions of slow rotation and weak magnetic field, the star is well approximated 
as a sphere. They found that toroidal fields are generated by magnetic braking, and the toroidal flelds and angular 
velocities oscillate independently along each poloidal field line. The incoherent oscillations on different fleld lines stir 
up turbulent-like motion on tens of Alfven timescales, a phenomenon called phase mixing (see |l2j[ and references 
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therein). In the presence of viscosity, the stars eventually are driven to uniform rotation, with the energy contained 
in the initial differential rotation going into heat. 

Most recently, Duez et al. 0, [3| and Shibata et al. performed GRMHD simulations on rapidly, differentially 
rotating magnetized neutron stars in axisymmetry using two, newly developed GRMHD codes 16l Il7l| . They found 
that if a star is hypermassive (i.e. the star's mass is larger than the maximum mass of a uniformly rotating star), 
magnetic braking and MRI will eventually induce collapse to a rotating black hole surrounded by a hot, massive torus 
with a coUimated magnetic field aligned along the spin axis. This system provides a promising central engine for 
short-duration gamma ray bursts T5]. They found the behavior of nonhypermassive, differentially rotating neutron 
stars to be quite different. If a star initially spins at a rate exceeding the limit for uniform rotation ( "ultraspinning" 
case), then instead of collapsing, such a star settles to an equilibrium state consisting of a nearly uniformly rotating 
core surrounded by a differentially rotating torus. Although this torus maintains differential rotation, the angular 
velocity is constant along the magnetic field lines, so further magnetic braking will not occur. However, if a magnetized 
star initially exhibits rapid, differential rotation at a spin below the limit for uniform rotation ("normal" case), the 
star settles into a uniformly rotating configuration. 

The purpose of this paper is to study the same magnetic effects examined by Duez et al. [T^L IT^ , but in the slow 
rotation, weak magnetic field limit (i.e., ^ T <C \W\). Analyzing the behavior of stars with such weak, but 
astrophysically realistic, magnetic fields on the magnetic braking (Alfven) timescale requires simulations spanning 
~ lO^M for the models we consider. Thus the primary challenge of this work is its exorbitant computational expense. 
To overcome this difficulty, we adopt a perturbative metric approach similar to the one developed by Hartle 18, 19], 
valid to first order in the angular velocity fl. The second computational challenge is to resolve the wavelength of 
the dominant MRI mode, which is small for the weak initial fields we wish to treat (Amri oc B). Solving the metric 
equations via perturbation theory allows us to adopt sufhciently high spatial resolution to track MRI. 

Two aspects of our perturbative approach make simulations significantly less costly than a nonperturbative, metric 
evolution at a given resolution and ultimately allow us to perform simulations at roughly 1 /4 the total computational 
cost. First, the perturbed metric is time independent except for the (/^-component of the shift, P'^, which gives rise 
to frame-dragging. The shift P'^ varies with f2 on the Alfven timescale, equivalent to many thousands of (Courant) 
timesteps in a typical simulation. Our perturbative metric solver uses a simple, ordinary differential equation (ODE) 
solver to compute the shift and allows us to skip many timesteps between matter evolution updates. Second, nonper- 
turbative metric schemes incorporate approximate asymptotic outer boundary conditions, which cause problems if the 
outer boundary is moved too close to the star. The perturbed metric on the other hand depends only on quantities 
defined within the (spherical) star, so the outer boundary of the MHD evolution grid may be moved significantly 
closer to the stellar radius with significant reduction in computational expense. 

The validity of our slow-rotation perturbative code is tested to ~ 2 Alfven timescales (^ lO'^M). We evolve a star 
with both the perturbative and the nonperturbative Baumgartc-Shapiro-Shibata-Nakamura (BSSN) pol | gravitational 
field scheme, and compare the results. We validate our simulations self-consistently by checking that the evolution 
data satisfy the slow rotation, weak magnetic field assumptions. 

We study resolution-dependent MRI effects using two techniques. First, we vary the grid spacing A at fixed initial 
magnetic field strength, and second, we vary the initial magnetic field strength (oc Amr t ) at fixed A. We observe 
MRI-driven rapid growth of the poloidal fields when Amri/ A ^ 10 (consistent with jlJillJI) and convergence in the 
growth rate if Amri/ A ^ 25. However, convergence in the maximum amplitude of these fields is not achieved even 
when Amri is resolved to ~ 41 points. This is due to small-scale turbulence intrinsic to MRI and to axisymmetry. 

In addition to an analysis of MRI, we also study magnetic braking. Our simulations indicate that the winding of 
magnetic fields due to magnetic braking saps a considerable fraction of the energy associated with differential rotation 
in roughly one Alfven timscale Ia^ regardless of resolution or metric evolution technique. Once the rotation profile 
becomes more uniform, the magnetic fields begin to unwind, pumping differential rotation energy back into the star. 

In the above simulations, we choose a differentially rotating star with an angular velocity profile that initially 
decreases with increasing distance from the rotation axis. In our final simulation, we evolve a differentially rotating 
star possessing an angular velocity profile that initially increases with distance from the rotation axis. Magnetic 
braking should occur both models, but MRI should not occur in the later case Our code yields the expected 

result. 

The remainder of this paper is structured as follows. In Section II, magnetic winding and MRI are explained 
qualitatively. Section III presents the mathematical and numerical framework for our simulations. Section IV outlines 
our initial data and numerical input parameters. Section V analyzes the validity of our perturbative metric solver. 
Section VI discusses our simulation results, and Section VII summarizes our conclusions. We adopt geometrized units 
in which G — c — 1. 
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II. QUALITATIVE OVERVIEW 
A. Magnetic braking 

In an infinitely conducting plasma (MHD limit), magnetic field lines are "frozen-in" to the fluid elements they 
connect. We evolve differentially rotating stars in this limit with initially purelypoloidal magnetic fields. Differential 
rotation causes the magnetic fields to wind toroidally on the Alfven timescale U|2^ 

R f B \-' f R \-'/'/ M 



where B is the magnetic field strength, va ~ B/y/Anp ~ B^/E?/iM is the Alfven speed, and R and M are the radius 
and mass of the star, respectively. As shown in |l3l |. in the early stage of the magnetic braking, the toroidal magnetic 
field -B"^ increases linearly with time according to 

B'^{t;w,z) = wB'''^tmB'{0;w,z)din{0;w,z) (i^w,z), (2) 



where = v"^ is the angular velocity, -y' = u'/m* is the matter three- velocity, and zu = y^x^ + y'^ is the cylindrical 
radius. The coordinates are set up so that the rotation axis is along the z-direction. 

The increase of B^ with time adds energy to the magnetic fields and saps the energy available in differential rotation 
Tbr, conserving total angular momentum. After roughly one Alfven time Ia, Idr is exhausted llOl ITll| and i?'^ 
reaches a maximum, so eventually the growth of B^ must deviate from the linear relation Although Ia depends 
on the initial magnetic field amplitude, the amplitude at which the field saturates does not [g. 

B. MRI 

MRI-induced turbulence occurs in a weakly-magnetized, gravitating body if the angular velocity decreases with 
increasing distance from the rotation axis 0, l2lj| . According to a local, linearized perturbation analysis in the 
Newtonian limit 0, 0, |^, this instability causes the poloidal field magnitude to increase exponentially with an 
e-folding time tmri independent of the seed field strength before saturating: 

TMRI = -2 TT, . (3) 



The wavelength of the fastest-growing mode, Amri, is given by 



where 



aw 



is the epicyclic frequency of Newtonian theory. For a star with 5 In 17/9 In ti7 ^ —1, we have tmri ^ 2/r2 ~ Pc, where 
Pc is the central rotation timescale, and Amri ~ 27rw^/r2. 

For configurations we consider, Amri/-R ~ 1/6, where R is the stellar radius. Therefore resolution on the order of 
R/A ~ 20 (where A is the grid spacing) is necessary to resolve the fastest growing MRI wave mode Amri- However, the 
onset of MRI results in a buildup of small-scale MHD turbulence, so the actual resolution requirements for fine-scale 
MRI modelling are much higher. 

C. Upper Bound on Magnetic Energy 

The magnetic energy A4 in a magnetized star undergoing magnetic braking increases as differential rotation is 
destroyed. In the slow rotation limit, the gravitational potential energy W and internal energy U of the star do not 
change significantly, so the maximum possible value of M for such a star, A^max, is given by 

A^max - Tor, (6) 
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where Tdr is the kinetic energy associated with differential rotation. Note that Tdr may be estimated at t = by 
constructing a rigidly rotating star with the same total angular momentum J as the differentially rotating star and 
computing the difference in kinetic energy between the two configurations. Eq. therefore provides us a way to 
estimate the maximum allowed magnetic energy a priori. 



III. BASIC EQUATIONS AND NUMERICAL TECHNIQUE 



In this section, we describe two methods to evolve the metric: the perturbative approach and the nonperturbative 
BSSN scheme. The perturbative approach fSec. lIIlXI) takes advantage of the fact that the system is nearly spherically 
symmetric. With this scheme, the evolution of the metric can be simplified considerably. The nonperturbative metric 
evolution approach fSec. IIIIB)| is the same as that used in The Maxwell and GRMHD equations are discussed 
in Sec. IIII 01 They are evolved with the same high-resolution shock-capturing technique as in [la |. 



A. Perturbative metric evolution scheme 



For a slowly rotating, quasi-stationary axisymmetric star, the rest-mass density po and pressure P differ from those 
of a spherical star to second order in rotation frequency Further, if the stress-energy tensor satisfies the circularity 
conditions (see Eq. (|B3|) and [33), we can choose a coordinate system so that the only off-diagonal component of 
the metric is the frame-dragging term gt^. In this approximation, the line element may be written to first order in 
rotation frequency Q and magnetic field strength |_B| as 

ds^ = -e'^'^^'^dt^ + e^^'^^dr'^ - 2uj{t, r, e)r'^ sin^ 6dtd(j} + r'^{sm^ dd^)^ + dO'^), (7) 

where r is the areal radius. Thus full determination of the metric requires expressions for e'^W, e^W, and uj{t,r,6) = 
—(3'^{t,r,6). The first two quantities comprise the time-independent components of metric, computed once for all 
time using the initial spherical P and po- However, w is a dynamical quantity that depends on the rotation profile of 
the star. It must therefore be recomputed as the star evolves in a quasi-stationary fashion. 

As stated before, the metric Q is valid only if the stress-energy tensor satisfies the circularity conditions. To simplify 
our calculation, we also require that the azimuthal momentum energy density associated with the electromagnetic 
field be small compared with those associated with the fluid. These conditions are satisfied if (I) the meridional 
components of the fluid's velocity are much smaller than the rotational velocity, and (2) the energy density of the 
poloidal magnetic fields is much smaller than the energy density of the fluid. Table summarizes these conditions, 
which are derived in Appendix IbI 



1. Computing the time-independent metric components 

For small fl, the equilibrium star is spherical (the deviation from sphericity is of order fi^), and the metric compo- 
nents e'^^''^ and e^'^''-' are independent of time (to order fi). They can be computed by solving the Oppenheimer-Volkoff 
(OV) equations 

^ = 47r.V(0 (8) 
dr 

dP{r) _ [p{r)+P{r)][m(r)^'Aiir^P{r)] 
dr r[r — 2m(r)] ' 

= [f - 2m(r)/r]-\ (10) 



with boundary conditions 



dv{r) _ 2[m(r) -f 47rr3p(r)] 
dr r[r — 2m{r)\ 



(11) 



po(0) = Pc = constant , (12) 
m(0) = , (13) 
lim v{r) = . (14) 
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We close the above set of equations via a polytropic equation of state: 

P = Kpl, r = 1 + l/n, (15) 

where K is the polytropic constant, F is the adiabatic index, and n is the polytropic index. Note that p is related 
to by p = po(\ + e) where e = -P/[(r — l)/0o] is the specific internal energy. In our perturbative scheme, v{r) and 
m(r) are frozen to their initial values and are not evolved with time. Note that outside the star, the diagonal metric 
components describe the Schwarzschild spacetime, with mass m{r > R) — M. 



2. The time- dependent shift term u)[t,r,6) 

What remains is to compute time-dependent quantity ijj{t, r, 9) = —(3'^{t, r, 6). Given the slow rotation assumptions 
summarized in Tabled the momentum constraint equation in the Arnowitt-Deser-Misner (ADM) formalism (Eq. 24 
in |23|) yields the following partial differential equation for uj{t,r,9) (see Appendix IXI for further details): 



^4 g^f. 



du! 
dr 



d 



r2 .\r.^Ode 



sm 



sm 



r r 



where j(r) = exp{— [A(r) + i/(r)]/2} and j'{r) = dj{r)/dr. 

Following we solve Eq. by expanding ui and Q in terms of associated Legendre polynomials: 



(16) 



u;{t,r,e) = ^P/(cos0)c^z(i,r), 
1=1 

oo 

n{t,r,e) = ^P/(cos0)f7Kt,r). 



(17) 
(18) 



Due to the assumption of equatorial symmetry (in addition to axisymmetry) , all even terms in the above expansions 
vanish. Substituting Eqs. ((T7|l and ((TH|l into Eq. l(TB|). we obtain the same radial equation for each I as Hartle (Eq. (30) 
of [H): 



1 d 

dr 



duJi 
dr 



Adj_ 
r dr 



.)/2^a + i)-2' 



^d3 
r dr 



(19) 



where fti{r > R) = 0. Our analysis of this equation in the limits r ^ and r ^ oo (see Appendix yields the 
following boundary conditions 



UJi{t,r)\r^O 



ni{t,0)+Ai{t), if? = l 

- ^[4p(0) + 3P(0)]r!3(t, Oy Inr, if Z = 3 
167r[4p(0) + 3P(0)]17/(t,0)r2 



Ci(t)r-'-\ 



3[;(? + i)-i2] 



otherwise. 



(20) 



(21) 



where Ci (t) and Ai (t) are determined (using the shooting method) at a given time t by matching the interior (r < R) 
and exterior {r > R) solutions at the stellar surface r = R±: 



dr 



uji{t, R+ 



oJi{t,R-) 
d 
dr 



(22) 
(23) 



For the models we consider in this paper, we find that contributions from modes above / = 5 are negligible, so we 
only calculate modes up to and including 1 = 5. 
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TABLE I: Assumptions made in the slow-rotation approximation (see Appendix |B| for a derivation). 



Orthonormal Component 


Velocity [max. average]* 


Magnetic Field [max. average]* 




v'^ — fir sin 9 


[^"^^^ < 1 [6 X lO-"] 


e 




yO 
Q.r 


< 1 [0.03] 


P'^' « 1 [4 X 10-^] 


f 






< 1 [0.04] 


[^^K « 1 [3 X 10"'] 
47rpo^ 



* The above nondimensional ratios are local in space and in time, so we compute a 
mass-density weighted average of these quantities at various times. We denote the 
maximum value (in time) observed in our simulations "max. average" . 



B. BSSN metric evolution scheme 



The line element for a generic spacetime is written in the standard 3+1 form as follows: 

ds^ = -a^df + 7y {dx' + 13' dt) (dx^ + p^dt), (24) 

where a is the lapse, /3* is the shift and jij is the three-dimensional spatial metric. We evolve the metric and the 
extrinsic curvature Kij using the BSSN formalism The BSSN evolution variables are: 






= ^ln[det(7,j)] , 


(25) 


% 


= e-^H, , 


(26) 


K 




(27) 






(28) 




= ■ 


(29) 



The equations for evolving these variables are given in pOl ]. For the gauges, we use the hyperbolic driver conditions ( |25l 
|2^) to evolve the lapse and shift. 

We adopt the Cartoon method [23| to impose axisymmetry and use a Cartesian grid. In this scheme, the coordinate 
X is identified with the cylindrical radius nj, the y-direction corresponds to the azimuthal direction, and z lies along 
the rotation axis. For example, for any 3- vector V, = V^, and = wV^ . 



C. Maxwell and MHD Equations 



In terms of the Faraday tensor F^^'^ , the MHD condition is given by 

F^^'u, = = 0, (30) 

where E'^^-^ is the electric field measured by an observer comoving with the fiuid. As in [l^. we evolve the following 
set of variables: 

P* = a77po"", (31) 

f = a'^T°°^p, , (32) 

= ay^T\ , (33) 

B' ^ VlB\ (34) 

where 7 = det(7ij ), Bf" = Fo^pn^, denotes the magnetic field measured by a normal observer, and is the unit 
normal vector orthogonal to the time slice. These variables satisfy the following evolution equations: 

9tU + V • F = S, where (35) 
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dtV - dt 



V-F = a 



f 



and 



S = 







The stress-energy tensor T^^'^ for a magnetized, infinitely conducting, perfect fluid is given by 



(36) 



(37) 



(38) 



(39) 



Here, h = 1 + e + P/po is the specific enthalpy, and y/ATrb^ — -^(u) magnetic field measured by an observer 

comoving with the fluid, which is related to i?^ by 



(«)' 



(40) 



where P^"' = g^"" + u^'u''. 

We evolve Eq. 1351) using a high-resolution shock-capturing scheme as in Specifically, we use the piecewise 

parabolic method (PPM) '28'| algorithm for data reconstruction and the Harten-Lax-Van Leer (HLL) flux formula 
for the approximate Riemann solver. 



D. Diagnostics 

During the simulations, we monitor the following conserved quantities: rest mass Mq, angular momentum J. We 
also monitor the ADM mass M , which is nearly conserved, as the energy emitted as gravitational radiation is negligible. 
We also compute the rotational kinetic energy T, magnetic energy M., internal energy t/, and gravitational potential 
energy W . All of these global quantities are calculated using the formulae given in |13l |. 



IV. INITIAL DATA AND NUMERICAL PARAMETERS 



To understand the behavior of slowly-rotating, weakly magnetized neutron stars, we perform four studies. First, in 
our "MRl Resolution Study," we start with a differentially rotating, poloidally magnetized configuration in which the 
angular velocity decreases away from the rotation axis. We then evolve this star at various resolutions, with the goal 
of uncovering the detailed, resolution-dependent behavior of MRL In our second study, the "P Variation Study," we 
evolve the same star as in the first study at lowest resolution, varying only the strength of the initial poloidal fields. 
This study also examines the resolution-dependent nature of the observed MRI by varying B and hence Amri at fixed 
spatial resolution. Finally, in the "Rotation Profile Study," we evolve the same star as with our "MRI Resolution 
Study," changing the angular velocity distribution so that it initially increases with distance from the rotation axis. 
In this study, we expect to observe magnetic winding, but not MRI (out to ^ Ua)- As a code test, we also perform 
the "Rigid Rotation Profile Study," where we explore the same configuration as with the first study, only with solid 
body rotation at the same total angular momentum J. We expect that the magnetic field will not change in time 
and have no effect on the star. Tables ITU and ITTll present a summary of initial parameters for the stars we consider in 
these studies. 

For simulations using the BSSN metric solver, we construct initial data for a differentially rotating, relativistic star 
in equilibrium using the code of Cook et al. |30j | with the following rotation law: 



(41) 
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TABLE 11: Initial models: Magnetic field-related parameters 



Study 




{tA)yM 


(tA)/ 


( M \ 




{tA)/Pc 


(|B|ma°x)7 


( l-4Mo ^ 




{M/\w\y 


r/|vi/|" 


\ 1.4MQ ) 




Rigid Rotation Profile 


6.1 X 10"'' 


4800 


33ms 


10.2 


4.9 X 10'*G 


4.4 X 10"'' 


4.55 X lO-" 


MRI Resolution 


6.1 X 10"'' 


4800 


33ms 


17.9 


4.9 X 10'*G 


4.4 X 10"" 


4.85 X lO"-" 


B Variation 


4.97 X 10"" 
1.96 X 10"^ 
6.1 X 10"^ 


16200 
8100 
4800 


112ms 
56ms 
33ms 


61.1 
30.7 
17.9 


1.4 X 10'*G 

2.8 X lO^'^G 

4.9 X 10^*G 


3.8 X 10"' 
1.5 X 10"® 
4.4 X 10"® 


4.85 X 10"'' 
4.85 X 10"^ 
4.85 X 10"^ 


Rotation Profile 


6.1 X lO"'' 


4800 


33ms 


7.2 


4.9 X 10'*G 


4.4 X 10"" 


4.90 X 10"'' 



C is the maximum value of at i = 0. 



' (tA) is the mass density- weighted Alfven time, given by Eq. 
* |Sj^x is the maximum magnitude of the magnetic field at t = 0. 
** A1/|VK| is the initial ratio of magnetic energy to gravitational potential energy. 
r/|VF| is the initial ratio of kinetic energy to gravitational potential energy. 



where R is the equatorial coordinate radius, f2c is the central angular velocity, and A is a constant parameter which 
determines the degree of differential rotation. In the Newtonian limit, this rotation law reduces to the so-called 
"j-constant" law: 

n = ^ . (42) 

1 + " 

For a slowly rotating star, the spatial metric is nearly conformally flat 7^ w ffj^^- The Cook et al. code uses 
spherical isotropic coordinates, so to obtain the desired 7ij, we only need to transform the Cook et al. initial data to 
Cartesian coordinates. 

For simulations with the perturbative metric solver, we set up the initial data by first computing the diagonal 
components of the metric and the hydrodynamic quantities by solving the OV equations lP^- (|ll|l . Then the shift 
is computed via the perturbative technique described in Section flll Al with the angular velocity distribution computed 
by either solving Eq. (|41|l in the slow rotation limit (the "Rotation Profile Study"), or using the solution of Eq. I|41|) 
as computed by the Cook et al. code ("MRI Resolution Study" and "i? Variation Study"). Note that since the 
initial data computed by the perturbative technique is only accurate to order the resulting star will undergo small 
amplitude oscillations [due to 0{VP') effects]. Further, to more easily compare perturbative simulation results with 
those using the BSSN scheme, we perform the coordinate transformations necessary to facilitate evolution of the 
Maxwell and MHD equations in the same Cartesian coordinates as in the BSSN evolution scheme. 

In our "MRI Resolution" and "S Variation" studies, we consider an A — \ differentially rotating star which satisfies 
the n = \ polytropic equation of state (EOS). Other parameters are set so that the equilibrium star possesses the 
following properties: the ratio of equatorial to polar radii Rp/R — 0.98, central rotation period Pc — 2Tr /flc = 264. 7M, 
compactness M/R = 0.182, ratio of angular velocity at the equator to that at the center fieq/^ic ~ 0.3, and = 
4.88 X 10"'^. The mass of this star is determined by the polytropic constant K, which we set to unity. However, our 
results can be easily rescaled to any values of K (see 113), and hence to any values of the mass. For example, the model 
we just described has R = 9.2 (M/1.4Mo)km, pc = 1.54 x lO^^ (1.4Mo/M)2g/cm3, and = 1.8 (M/1.4M0)ms. 

The "Rotation Profile Study" involves the same star as in the "MRI Resolution Study" , but with rotation profile 
parameters set so that A = —1, which corresponds to il.{R)/il.c — 2.8. 

Next, we add a small seed magnetic field to the stellar models above by specifying the vector potential Ai = A^/iSi'^ 

as 

A^ = max[Ab{P - Pcut), 0] , (43) 

where the pressure cutoff Pcut is set to 4% of the maximum pressure (Pcut = 0.04Pmax)- The strength of the initial 
seed field is determined by the constant Ah and may be characterized by the parameter C, the maximum value of 
6VP at t = 0. 

The strength of the magnetic field can also be measured by the mass density-averaged Alfven time (i^) defined as 



R 



1 f .3 r' 

VAP*d X 



RMo 



(44) 



where va = yW/ifioh + W) is the Alfven speed. Since Amri is a local quantity, we define the magnetic energy 
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TABLE III: Initial models: Parameters related to resolution and numerical evolution 



btuay 


Metnoa 


(K/A) 


o 




(tstop /-rc) 


Rigid Rotation Pronle 


Perturbative 


100 


6.1 X lO"'' 




10.2 (ns) 




Porturbative 


75 


6.1 X 10^*" 


12.2 


35.8 (ns) 






100 




16.3 


15.5 


MRI Resolution 




150 




24.4 


13.1 






OAA 
ZUU 




oZ.D 


io.o 






zoU 




4U. ( 


iz.o 




Nonperturbative 


75 


6.1 X 10 


12.2 
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30.7 (ns) 


B Variation 
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Nonperturbative 


75 
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6.1 X 10"^ 


12.2 
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Rotation Profile 
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* f? is the equatorial coordinate radius of the star, and A is the grid spacing. 
^ tstop is the time at which the simulation was stopped due to loss of accuracy, 
which happens soon after the magnetic field hits the outer boundary, (ns) 
indicates that the magnetic fields have not hit the outer boundary at 2{tA), and 
the simulation was terminated at the indicated time. 



density-weighted average of Amri as follows: 



^^^^^^ = Jb^^d^x ' 



where 



r r \ j ^Mm{r) if < rMRi(r) < tA , 
AmriW = <„ ^, . • (46) 

I otherwise 

Here Amri and tmri are calculated by Eqs. and Q, respectively. The cutoff in Eq. H46() is set so that we only 
consider the region where the MRI is present (tmri > 0) and where the MRI timescale is less than the Alfven time 
(where tmri > ^a, magnetic braking is expected to dominate). 

All our nonperturbative simulations are performed on a square grid with outer boundary at 2.0R (= 11 Af). Our 
perturbative simulations on the other hand use an outer boundary of 1.2R, with metric updates every 8-10 timesteps. 
We have verified that if the outer boundary is set to 1.5R instead, all quantities we studied are the same to within 
~ 1% until the magnetic field hits the 1.2R outer boundary. Due to this loss of accuracy, we stop our simulations 
soon after this boundary crossing. The time at which each simulation was stopped, <stop, is listed in Table HTll By 
~ 2tA, both magnetic winding and MRI are fully developed. 

We specify resolution by the quantity R/A, where R is the stellar radius and A is the grid spacing. Thus a 
simulation with 150^ points and outer boundary at 2R has R/A — 75, and a simulation with 90^ points and outer 
boundary at \.2R has Rj A = 75 as well. In these simulations, MRI does not become evident until t ~ 6Pc (e.g., 
see Figure EJ. Thus for computational efficiency in our highest resolution run, R/A = 250, we evolve the star at 
resolution R/A = 100 until t = bPc and then rcgrid to R/A = 250. Table Hill summarizes the resolutions chosen in 
our simulations. 



V. CODE TESTS 



A. Test of the perturbative shift solver 



To verify that our perturbative shift solver produces the shift P'^ accurately to order il, we compute (3'^ for 
differentially rotating, equilibrium star models and compare the results with those computed without approxima- 
tion by the Cook et al. code [s^- Figure ^ shows the error, (5/3"^, for three models with the same central density 
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FIG. 1: Relative error in /J*^ along equatorial plane of star at t = 0. Exact solutions from Cook et al. code [3Cj are compared to 
perturbative shift solver results at T/\W\ = 2.43 x 10"^ (solid black lines), 4.9 x 10"^ (dashed blue lines), and 9.9 x 10"^ (dotted 
red lines). Here 5/3*^ = {P^^^i^ — /3pert)/[(/^Cook + /3pcrt)/2]- To demonstrate the approximate scaling 5/3*^ oc r/|W|, we have 
14- c fcr-fr^r r^f A w f>,o ^»co T/lwl - 2.43 X 10^=^, aud bv a factor of 2 for the case T/\W\ = 4.9 x 10"^. 



multiplied the value of 5/9'^ a factor of 4 for the case r/| W| 



L = 



0.0 (g 





FIG. 2; Snapshots of the poloidal magnetic field lines (contours of A^) at t/ {Ia) = 0, 0.5, and 1.0 for the uniformly rotating 
star. The field lines are drawn for = ^0,min + (^0,max — ^c/),min)*/20, (i = 1-19), where ^0,max and A^^^^in are the maximum 
and minimum values of A^, respectively, at the given time. The dashed line in each plot indicates the initial stellar surface. 



(pc = 1.54 x 10^5 (1.4M0/M)2 g/cm3) but with various T/\W\. We see that SP'I' decreases as fi^ cx T/\W\, 
expected. This shows that our perturbative shift solver accurately calculates P'^ to first order in fl. 



B. Rigid Rotation Profile Study 



When the star is uniformly rotating, a (weak) poloidal magnetic field should not change with time and it should 
have no effect on the star. To test our code, we evolve a uniformly rotating star with the physical parameters specified 
in Tables HTl and UTTl ( Rigid Rotation Profile Study). We follow the star for one Alfven time ((t^) = 10.2Fc = 4800Af). 
Figure [3 displays snapshots of the poloidal magnetic field in time, and Figure O shows the evolution of the rotational 
profile on the equatorial plane. As expected, neither the star's rotation profile nor its magnetic field change significantly 
over the Alfven timescale. 
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FIG. 3: Snapshots of the rotation profile at the equatorial plane at t/{tA) = and 1.0 for the uniformly rotating star. 




FIG. 4: |-B^|max vs. time, with perturbative results plotted on the left and nonperturbative (BSSN) on the right. Resolutions 
of R/A = 75, 100, 150, 200, and 250 are shown with solid, dotted, short dashed, long dashed, and long dash-dotted lines, 
respectively. The short dash-dotted line represents an approximate slope 7 — 1/(5. OPc) for the exponential growth rate of the 
MRI, SB'' oc e^*. Note that {tA)/Pc = 17.9. 



VI. NUMERICAL RESULTS 



A. MRI Resolution Study 



In the MRI Resolution Study, we perform simulations on a magnetized differentially rotating star (see Table ^Oj 
at resolutions R/A = 75, 100, 150, 200, and 250. The initial magnetic field is set so that C = 6.1 x 10"^. To 
demonstrate that the two metric solvers yield the same results, we also perform simulations with the BSSN metric 
solver at the two lowest resolutions (-R/A =75 and 100). Figure 01 displays the maximum magnitude of \B^\ as 
a function of time. MRI causes the sudden increase of l-B'^lmax at i w dPc- As resolution is increased, jiJ'^lmax 
saturates at larger values. Due to the turbulent nature of the MRI, we do not achieve convergence, even at the 
highest resolutions (-R/A = 200 and 250). However the exponential growth time of the MRI, tmri, does converge 
at the highest resolutions (-R/A = 150 250). The numerically determined value for tmri is ~ 5. OPc: which does 
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FIG. 5: Evolution of |-B''|max(t). Simulations with resolution i?/A = 75, 100, 150, and 250 are shown with solid red, dotted 
blue, short dashed cyan, and long dashed green lines respectively, with perturbative results on the left and nonperturbative 
(BSSN) on the right. The dash-dotted black line represents the expected early-time linear growth of |-B^|max, as predicted by 
Eq. 0. 



not significantly deviate from the linearized, Newtonian theory estimate (obtained by applying Eq. P|) at t = 0) of 
TMRi.min — 5.7Pc- We See that regardless of resolution or metric evolution scheme, the MRI-induced amplification of 
the magnetic field above its initial value becomes evident hy t fa 6Pc- 

FigureElplots the maximum value of as a function of time. The straight lines in each plot indicate the expected 
growth rate via magnetic braking in the hnear regime, according to Eq. 10). Notice that at early time, |i3^|max agrees 
well with the expected linear growth. However, the slope begins to flatten once the magnetic field becomes strong 
enough to induce fluid back-reaction. Later, the slope of |i?^|max increases once again before flattening at the point 
of saturation. As resolution is improved, the saturation point in |i?^|max occurs at earlier time. We find that the 
sudden increase in the slope correlates well with the increase of |-B^|max in the MRI plot of Figure 0| We also find 
that at the same time, the region at which |B^|max occurs shifts to the region where the MRI occurs. This is further 
evidence that the sudden increase in |i?^|max results from MRI-induced magnetic field rearrangement. 

To further explore magnetic field rearrangement in our star at different resolutions. Figure |S1 provides snapshots 
of poloidal magnetic field lines (i.e., contours of the vector potential A^) at various times. Notice that MRI induces 
the largest distortion of the field lines in the outer equatorial region of the star. This is consistent with the linear 
analysis: Eq. Q, together with the star's angular velocity profile, gives a shorter tmri near the outer part of the star. 
Similar behavior has been observed in simulations of magnetized, rapidly and differentially rotating stars [Tslfl^ . The 
distortion becomes more prominent at finer resolution, which suggests that more and more small-scale MRI modes 
are being resolved as resolution is improved. 

Next we analyze how the rotation profile of the star changes as a result of its magnetic field. Figure [7| shows 
the equatorial rotation profile at various times. Consistent with the results of Q, we find that magnetic winding 
destroys the differential rotation profile on the Alfven timescale, causing the star to rotate nearly as a solid body with 
n = f^const at t « l(iA)- Here ficonst is the angular velocity of a uniformly rotating star with the same rest mass 
and angular momentum as the star under study. When the toroidal field saturates around t « l{tA), the magnetic 
fields begin to unwind, eventually causing the rotation profile to increase with increasing radius at roughly 2{tA)- In 
addition, the MRI stirs up a turbulent-like flow, causing the bumpy rotation proflle seen at later times. 

As the magnetic flelds are wound, the magnetic energy M saps kinetic energy T associated with differential rotation 
(Fig. IS)) until the star rotates nearly as a solid body. At this point, we find that the star's kinetic energy sinks to 
its minimum value shortly after Ai reaches maximum. Note that, consistent with Fig. 0] the maximum of Ai occurs 
earlier and earlier as resolution is improved due to an interplay between MRI and magnetic braking. After M reaches 
maximum, the fields unwind, pumping energy back into differential rotat ion, as shown at lowest resolutions in Fig. |S1 
We speculate, based on the a-disk model jU and on our previous work that the oscillations of T and M 

will continue for many Alfven times until the rotational kinetic energy associated with differential rotation is dissipated 
by phase mixing caused by MRI-induced turbulence. However, since the star is slowly rotating (T/|VF| = 4.88 x 10^^) 
with weak magnetic fields {Ai <C T), the Alfven time is long, (t^) = 10.2Pc = 4800M. It is therefore computationally 
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FIG. 6: Snapshots of poloidal magnetic field lines at various times, with resolution R//S. — 75 (top), 150, (middle), and 250 
(bottom). The field lines are drawn for = A^^^nin + {A^^max — A,/,,niin)i/20, (i = 1-19), where v4^,max and A^^min are the 
maximum and minimum values of A^, respectively, at the given time. The dashed line in each plot indicates the initial stellar 
surface. 




FIG. 7: Rotation profile fl measured in the equatorial plane at various times, comparing perturbative scheme (left) with 
nonperturbative (right) at R/A = 75 resolution. The straight line indicates the solid body angular frequency ficonst of a star 
with the same angular momentum J and rest mass AIo- 
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FIG. 8: Rotational kinetic (T) and magnetic (M) energies vs. time. For a given energy E, we define AE = [E{t) — E{0)]/T(0). 
The left plot compares results from the perturbative (solid red) and BSSN (dashed blue) metric algorithms at R/A = 75 
out to t = 2{tA), and the right plot shows the same at R/A = 100 (perturbative: dotted red, BSSN: dashed blue) and 250 
(perturbative only: solid green). The rotational kinetic energy of rigid body rotation for the given J is plotted at the bottom of 
each graph. The data for the perturbative runs have been smoothed to remove unphysical oscillations (see the text for details). 



taxing to accurately evolve the star for many Alfven times, even if the perturbative metric solver is used. 

As discussed in Section IIVI perturbative metric solver initial data are only accurate to order f2. This causes 
oscillations to arise in our perturbative T data at the level of AT « 0.0065, where AT = [T - T(0)]/T(0) is the 
fractional deviation of the rotational kinetic energy from its initial value. The oscillations are evident in raw Fig. |S| 
perturbative data. A simple, local (in time) averaging technique is used to smooth out these oscillations in Fig. |S1 
Although they are smaller by an order of magnitude, we remove the oscillations in our perturbative A4 data as well, 
for consistency. Notice also in Fig. |S1 that the value of A^max is comparable to that derived in Eq. ©. 

Figure|H|demonstrates that the angular momentum J in our long-term simulations is well-conserved out to 2{tA)- We 
see that angular momentum is lost at nearly a constant rate in our perturbative simulations, but the loss decreases with 
increasing resolution. In addition to angular momentum conservation, the binding energy Mq — Madm is conserved 
to within « 0.5% in perturbative and « 2.5% in BSSN simulations. 



B. B Variation Study 



In this study, the grid resolution is fixed at R/A = 75, and only the strength of the initial magnetic field is varied. 
We perform both perturbative and BSSN simulations. The MRI wavelengths in these simulation are (Amri)/A =3.6, 
7.2, and 12.2. The corresponding magnetic field strength parameters are given in Table Ull Note that the last case is 
the same as the lowest resolution case in the "MRI Resolution Study" . 

In Fig. ^1 we plot |i3^(i)|max for the three magnetic field strengths. Although (Amri) depends on initial magnetic 
field strength (Eq. l|^ ) , the MRI e- folding timescale tmri does not (see, e.g. Eq. ||3Jl). MRI is observed only when 
(Amri)/A > 12, which is consistent with the results of [islIT^. 



C. Rotation Profile Study 



In this study, we consider a differentially rotating star in which Q increases with cylindrical radius. We expect to 
see magnetic braking but not MRI in simulations of this star. 

Figure ITTI presents the same plots as those in Figs. 0HZI for this star. The top left plot indicates that, as expected, 
MRI is absent from this simulation. Although magnetic braking does appear (top right) as predicted, the curve is 
nearly parabolic and does not exhibit a sudden increase in slope as observed in Fig.jS] This further supports the notion 
that the sudden increase in slope in Fig. [S] results from an interplay between magnetic braking and MRI. Note also 
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FIG. 9: Relative change in angular momentum J (AJ = [J — J(0)]/J(0)) vs. time at resolutions i?/A — 75 (solid red), 150 
(dashed blue), and 250 (dotted green). Recall that the R/A — 250 data is from a regridding run, where the data before 
t/{tA) ~ 0.3 is at R/A — 100, and 250 after. This explains the varying slope in the R/A = 250 data. We plot results from the 
perturbative scheme only; the nonperturbative technique conserves J to the same degree. 
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FIG. 10: IB^Imax vs. time, with different initial magnetic field strengths, so that (Amri)/A = 4.0 (bottom line), 8.0 (middle 
line), and 12.2 (top line). Resolution is fixed at R/A — 75. Data from the perturbative spacetime evolution method are shown 
in the left plot, and BSSN method in the right. 



that the MRI- induced magnetic field distortion does not appear in the poloidal plane (middle three plots). Finally, we 
see that since the magnetic field-shifting effects of MRI are absent and the magnetic field is initially confined to the 
high density inner region of the star (Fig.|SJ), magnetic braking is incapable of fiattening the rotation profile (bottom 
plot) in the less dense outer layers of the star. 



VII. CONCLUSIONS 



Our perturbative metric evolution algorithm yields results quite similar to those produced by the nonperturbative, 
BSSN-based evolution scheme. Because the outer boundary may be moved inward and the metric updated on a 
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FIG. 11: Results from the Rotation Profile Study. The plots from top left to bottom are analogues to Figure 0J|7| Top left: 
|-B^|max vs. time, top right; |i3^|max vs. time, middle three plots: magnetic field lines (contours of A^), bottom: equatorial 
rotation profile. 



physically relevant timescale, perturbative metric simulations may be performed at ~ 1/4 the computational cost 
of BSSN metric simulations in axisymmetry for slowly rotating, weakly magnetized equilibrium stars. However, we 
found that loss of accuracy did occur after magnetic fields hit the outer boundary. This problem could be efficiently 
solved in future work by extending the outer boundary further from the star. 

The stars we study in this paper are weakly magnetized ^ 10~^) and slowly rotating (T/|VF| w 0.005). As 

a result, the Alfven timescale is very long (~ 5000Af), and accurate simulations spanning many Alfven times would be 
prohibitively expensive. In this paper, we reliably evolved such stars out to two Alfven times at moderate resolution 
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(i?/A = 75). In these simulations, we found that magnetic fields wind and unwind on the Alfven timescale, resulting 
in a trade-off between the kinetic energy in differential rotation and magnetic energy. Since we could not perform 
simulations spanning many Alfven timescales while accurately resolving MRI, we can only speculate, based on the 
a-disk model and our previous work that the oscillations between the magnetic and kinetic energy will 

be damped by MRI-induced dissipative processes over many Alfven times. 

In order to observe MRI in our simulations, we find that the spatial resolution must be set so that Amri/A > 10, 
in agreement with the results of 0,^3- have also verified that MRI is not present if the star's angular velocity 
initially increases with increasing distance from the rotation axis (i.e., the MRI is absent when d^fl > 0). 

We find that as resolution is increased, the effects of MRI become more and more prominent. Due to the turbulent 
nature of MRI, we do not achieve convergence of field amplitude, even if Amri/A = 32.6 or 40.7. However, we found 
the e-folding time of MRI, tmri does converge. The numerically determined value, tmri ~ S.OPc is consistent with 
the value predicted by the linearized, local Newtonian analysis (5.72Pc)- The small difference is due in part to the fact 
that our star is relativistic. In addition, the linearized analysis assumes that Amri is much smaller than the length 
scale on which the magnetic field changes (i.e. Amri ^ B/\WB\), which is not quite satisfied in our magnetic field 
configuration. 

Finally, we note that the behavior of the MRI is expected to be different in a full 3D calculation because of the 
effect of nonaxisymmetric MRI induced by a toroidal magnetic field. Turbulence may arise and persist more readily 
in 3D due to the lack of symmetry. More specifically, according to the axisymmetric anti-dynamo theorem p3ll l32j| , 
sustained growth of the magnetic field energy is not possible through axisymmetric turbulence. Thus proper treatment 
of MRI in differentially rotating neutron stars requires high resolution simulations performed in full 3-1-1 dimensions. 
The computational cost of such simulations with existing 3-1-1 metric evolution schemes has thus far been prohibitive, 
but with our perturbative metric solver it may be possible to perform 3-1-1 simulations of weakly but realistically 
magnetized, slowly rotating stars at a small fraction of the computational cost. 
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APPENDIX A: DERIVATION OF SHIFT EQUATION (EQ. 

Here we derive the equation for the shift (Ea. I16|l . starting from the usual momentum constraint equation in the 
3-1-1 ADM decomposition of Einstein's field equations. Recall from Eq. {Tj) the perturbative line element is given by 

ds^ = -e'^^^'^dt^ + e^^^'Ur^ + 2P'''{t,r,0y sin^ 0dtd<f) + r^isin^ ed(l)'^ + de^) + 0{n^) + 0{B^). (Al) 

We begin by writing the momentum constraint equation (Eq. 24 in 24]) 

DjK^, - D,K = 87rj„ (A2) 

where K — j and ji — —j^iUcT'^b- It can easily be shown from the ADM 3-metric evolution equation (Eq. 35 in 
H) dt"fij = —2e'^^^Kij + Cp^ij that in this approximation Kij is given by 

= ^a^J = -2e'''^K,, + + (A3) 

where (3i\j denotes the spatial covariant derivative and the first equality reflects the stationarity of the 3- 

metric in this approximation. Taking the trace of this equation and applying the identity y = -^[y^(3^) j (where 
7 = detjij) yields an expression for K: 

_2e-/2if+ A(^A-),^. =0. (A4) 

Note that (^/7/3-')j — {y/lP'^),<p = by axisymmetry, so Eq. HA2p may be rewritten 

Kh\j = STTj,. (A5) 
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It follows from the definition of ji and the metric ljAl|l that ji = e^/'^T'^i. We split the stress-energy tensor T^^, as 
well as ji, into a fluid part and an electromagnetic part: 



rp _ rpF I ^EM 



T. 



EM ^ ^2 

J 



F _ l//2rTnF0 

i — e 1 i , 



J 



EM _ ^ij/2jiF.M0 



(A6) 
(A7) 

(AS) 
(A9) 



Our approximation as summarized in Table and derived in Appendix ^ensures that is the dominant component 
of ji. In this approximation, — (0, 0, fl) and 



Thus i = 4> remains the only nonzero component of Eq. IjASp to first order in and B: 

K^^\j = ^-Kpohe-^l^in + P'>')r^ sin^ 6 . 

Next we expand our expression for K'^ 

= ^(V7r2sin2 0if^-^),,-r'=^„i^"'7feZ. 



(AlO) 
(All) 

(A12) 
(A13) 



It follows from Eq. and the identity = ^ + ^kfi^ that 

2e-/2x»^' = + r^fc/3'=)7J"' + + (A14) 
After computing the necessary Christoffel symbols, we find that the nonvanishing components of K'^^ are: 



1 



Plugging these expressions into Eqs. ljA13|) and then ljAll|) yields the equation governing the shift: 



■ e 2 



-[sin" 0LJ.g].g + -j\r)u; ^ -j'{r)n, 



sin^ y ' r r 



(A16) 



(A17) 



where lo 



and j{r) = e 



Finally, we perform the angular decompositions as in Eq. p7f) and Eq. p8fl and substitute them into Eq. IjAlTp to 
obtain the equation 



r 



(A18) 



This equation is the same as Eq. (30) in 19], which was derived using a different approach. Each w; must also satisfy 
the boundary conditions at the origin and at infinity (see below), as well as the matching conditions [Eqs. 1)22(1 and 
(j2SJ] at the surface of the star. 



1. Boundary condition at the origin 



We now determine the boundary condition at the origin by analyzing the terms of Eq. ((A 18(1 in the r ^ limit. 
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First we analyze j{r) and its derivatives, starting with j'(0) = — ^i(0)[A'(0) + J^'(O)]. From the OV equations, we 
obtain 



A'(0) 



^'(0) 



2[rm'{r) — 'm{r)] 



(r - 2m{r)y 
2[m(r) +47rr3p(r)] 







r[r — 2m(r)] 



(A19) 
(A20) 



where we have used the expression m(r) = (4/3)7rr'^p(0) as r 0. Thus j'(0) — since j{r) — e 2 is regular at 
r — 0. Next we examine j'{r)/r: 



/(r)/r|_o = /'(O) = {e-^r\r^o 

= -^j(0)(A"(r) + ^."(r))|,,^o 



1 , . d (2[rm'(r)-Tn(r)] 2[m(r) + Awr^ P(r)] 



dr 



[r — 2m(r)]^ 



2m(r)] 



r— >0 



- --^j(0)[4p(0) + 3P(0)]. 



(A21) 
(A22) 

(A23) 

(A24) 



We may therefore write Eq. (|A18|) near r ^ as follows 



r j(0) j(Oj 



(A25) 



with solution given by Eq. (|20|) . 



2. Boundary condition at infinity 

Outside the star, P — p — Q and the time independent (diagonal) metric becomes Schwarzschild, so = e^^ 
1 — 2M/r. The ODE governing the shift outside the star is therefore given by 



^il'^^-A^)]' - T^i^^^^^-^M = 0. (A26) 



Since 2M/r ^ 1 in the limit r — > 00, we may write Eq. HA18|I as 



1 d ^^dhJi^ l{l + l)-2 
dr dr 



^' = (A27) 



with solution given by Eq. 

Note that the analytic solution for Eq. I|A26|I exists and is given in terms of the hypergeometric function: 




if Z = 1 , and 

uJi = { (A28) 
;2Fi{l + 2,l- \;2l + 2\2M/r), otherwise. 



3. Rigid rotation case 



In the case of solid body rotation {Vl{r,6) = constant), only the I = 1 mode in Eg. 1191 contributes to H,. Thus the 
right-hand side of Eq. l|A18p is zero for I > 1. For / > 1, the solution uji — satisfies the boundary conditions at the 
origin and at infinity [Eqs. (|20|l and l|21() ] and the matching conditions at the stars's surface [Eqs. (|22|l and H23() ]. so 
uji — is the solution for / > 1. This coincides with the result cited in [T^ . 
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APPENDIX B: DERIVATION OF SLOW ROTATION APPROXIMATION INEQUALITIES 



In this section, we derive inequalities that must hold in order for our primary assumption in Appendix ^ [leading 
to Eq. l|ITO|) ] to be vaHd. 

We have assumed that the metric can be written in the form {T)) at all times, with the shift (3'^ = — w being the 
only non-diagonal component of the metric. For this to be true, the system has to be (approximately) stationary, 
axisymmetric and the stress-energy tensor has to be circular or nonconvective |33j| : 



(Bl) 
(B2) 



where e = d/dt and ^ — d/d4> are two Killing vector fields associated with stationarity and axisymmetry, respectively. 
These circularity conditions are satisfied if the momentum currents in the meridional planes are negligible compared 
with the axial component. Hence we require that 



Ij^I » \jf\ and Ij^I » , 

where the "hats" denote the orthonormal components. 

We split the stress-energy tensor into a fiuid part and an electromagnetic part: T^'^ — Tp" 



^EM' where 



2^ 



From this we obtain the following expressions for T^i^ and T^i^^: 



r 



nO F 



EM 



— pohu^Ui 

= b^u^u,-b^b, 



(B3) 

(B4) 
(B5) 



— [^,^B'-B.,B,) , 



47r 



(B6) 



(B7) 



where we have used the fact that + (3^ ~ and au^ = 1-1- O(ri^) for slowly rotating stars. Here the lapse a 

Since ji = e'^^^T^i, we may split the momentum current density in the same way: ji = jj + jf"^- Using the above 
expression for and the metric Q), we find 



3t 



F _ „v/2rp0 F 



■'4, 


- e-^/^pohrn 








< -'e 


\nr) 




< 







(B8) 



In the absence of magnetic fields, ji = and the conditions ljB3|l yield |f | ^ \fl\r and \v^\ <ti \^\r [Assuming 



-A/2 



0(1)]. Applying these inequalities, the electromagnetic part of ji becomes 



•EM 



Ji 



EM 








B^ 






B^B^ 


B^ / 


52 



(B9) 



For simplicity, we have ignored the magnetic field terms when computing the shift in Appendix ^ For this to be 
valid, we need to impose an additional condition: 
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Equations (jB8|) and ljB9|l together with the conditions ljB3|) and (|B10|I yield the following inequalities which must 
be satisfied for the shift equations in Appendix fXl to be valid: 
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Or 
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e\2 



< 1 



<C 1 , and 



< 1 



(Bll) 



For information on how well these inequalities are satisfied in our simulations, see Table I. 
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